Cells, cancer, and rare events: homeostatic metastability 
in stochastic non-linear dynamics models of skin cell proliferation 
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A recently proposed single progenitor cell model for skin cell proliferation [Clayton et al, Nature 
446, 185 (2007)] is extended to incorporate homeostasis as a fixed point of the dynamics. Unlimited 
cell proliferation in such a model can be viewed as a paradigm for the onset of cancer. A novel 
way in which this can arise is if the homeostatic fixed point becomes metastable, so that the cell 
populations can escape from the homeostatic basin of attraction by a large but rare stochastic 
fluctuation. Such an event can be viewed as the final step in a multi-stage model of carcinogenesis. 
This offers a possible explanation for the peculiar epidemiology of lung cancer in ex-smokers. 

PACS numbers: 87.17.Ee, 87.18. Tt, 87.19.xj 



The epidermis is the outermost part of the skin barrier. 
It comprises 10-20 layers of skin cells, which are predom- 
inantly keratinocytes fT] . Keratinocytes prohferate in the 
basal layer, move up through the middle layers, and are 
finally shed from the outermost layer at a desquamation 
rate of the order 10^ cells hr"'^ mm~^ 0. This process 
means that cells have to proliferate continuously in the 
basal layer to replenish the supra-basal layers. Faulty cell 
proliferation in the basal layer has important health-care 
consequences, for example basal cell carcinoma is one of 
the most prevalent forms of cancer 

Motivated by elegant in vivo experiments on labelled 
keratinocyte clones, Clayton et al. 0, d, @| recently pro- 
posed a novel single progenitor cell (SPC) model for basal 
layer proliferation. In the present paper I examine an ex- 
tension of this SPC model to include autoregulation, and 
homeostasis as a dynamical fixed point. Interestingly, 
something resembling the onset of cancer (i. e. carcino- 
genesis) arises naturally in the new model if the homeo- 
static fixed point loses stability in the direction of unlim- 
ited growth of the cell populations. One intriguing possi- 
bility is 'homeostatic metastability', in which a large but 
rare stochastic fluctuation results in the cell populations 
escaping from the homeostatic basin of attraction. Such 
a rare escape event can be viewed as the final step in a 
multi-stage model of carcinogenesis [tI, Q ■ 

In the SPC model, there are two cell types, or 'com- 
partments': progenitor cells A, and post-mitotic cells B. 
These proliferate according to 



A ^ A + A rate Ai, A 
A ^ B + B rate A3, B 



A + B rate A2, 
rate F. 



(1) 



The first three processes represent possible progenitor 
cell division pathways. The last process represents post- 
mitotic cells leaving the basal layer. Obviously, the pro- 
genitor cell division pathways must be finely balanced 
otherwise the cell populations either grow or vanish ex- 
ponentially [9]. Making the assumption therefore that 
Al = A3, Clayton et al. [3, HI write Ai = A3 = Ar and 
A2 = A(l — 2r). For mouse skin, they find A « 0.16 day^^ 



for the overall progenitor cell division rate, and 2r « 0.16 
for the branching ratio into the symmetric division path- 
ways (in later work 2r « 0.4 Q but the actual value is 
not too important). A steady state condition for the cell 
populations is 



Xp = F(l - p) 



(2) 



where p = rij^/n w 0.22 is the fraction of progenitor cells, 
n — UA + nB is the total cell density, and tia and tib are 
the individual cell densities. From this one determines 
F « 0.045 day~^, compatible with the above-mentioned 
desquamation rate. 

A phase space portrait and a representative stochastic 
trajectory for the SPC model are shown in Fig.[TJa). The 
trajectory was generated using a standard kinetic Monte- 
Carlo algorithm [lo'|, starting from a patch of area A 
comprising initially N ~ An = 200 cells, and assuming 
the cell types remain well-mixed. The actual value of A 
is irrelevant since the processes in Eq. ([T]) are all first 
order. In Fig. [1] the ordinate is n/no = N/Nq where uq 
(Nq) is the initial cell density (number). 

Whilst providing an eminently satisfactory explana- 
tion for the keratinocyte labelling data, the original SPC 
model includes a fine-tuning assumption, Ai = A3, mak- 
ing it structurally unstable in the language of dynam- 
ical systems theory ll]. As a consequence the model 
lacks homeostasis in the sense that it possesses a line 
of stable fixed points, shown in Fig. [TJa). Also it is 
not generally robust against perturbations, for exam- 
ple it cannot accommodate the introduction of a small 
population of stem cells jT2|, without making some ad- 
ditional fine-tuning assumptions. One obvious solution 
to this is to suppose that the cell division rates depend 
on the cell population densities ha and ns, representing 
the fact that cellular fates are governed by integration 
of autocrine and paracrine signalling factors 12, H, |lj]; 
indeed this idea was already suggested by Jones and Si- 
mon as an avenue for further investigation [l5| . In such 
an autoregulating SPC (ASPC) model, the fine-tuning 
would arise as a consequence of the cell population dy- 
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namics driving the system to a homeostatic fixed point. 
An ASPC model would be structurally stable from the 
point of view of dynamical systems theory, and also able 
to withstand perturbations, such as the presence of a 
small number of stem cells. 

To develop such an ASPC model, I start by introducing 
a control parameter q to describe a possible bias in the 
symmetric cell division fates, thus 

Ai-Ar(l-g), A2 = A(l-2r), A3 - Ar(l + g). (3) 

The parameters A, r and q replace Ai, A2 and A3 without 
loss of generality. The steady state conditions are given 
by g = and Eq. ([2]). From this point onwards, I delib- 
erately adopt a 'physics-oriented' approach in which the 
model is kept as simple as possible to expose the general 
mechanisms at work. In particular it is not claimed that 
the biology is accurately represented. The basic idea is to 
introduce a minimalist dependence of q and A on the cell 
population densities, to represent the integrated effect 
of the intercellular signals. I start by making q = q{p). 
Then the steady state condition selects the value of p for 
which q{p) = 0. I additionally suppose that A = A(n) is a 
decreasing function of the total number density (A' < 0) 
representing the fact that the progenitor cell proliferation 
rate should be reduced if the overall cell density increases 
[l^ . With these assumptions, fixed points of the dynam- 
ics are determined by q(p) = and A(n) = T[l — p)/ p. It 
is a straightforward exercise to show, in the language of 
dynamical systems theory, that a fixed point is a stable 
node if g' > 0, and a saddle if g' < 0. 
A concrete model of this type has 



A = An(^ 



g = tanh 



/3po(l - Po){p - Po) 



P(l - P) 



(4) 



where Aq = r(l — po)/po- This model has a stable node 
at a target population density n — no and progenitor cell 
fraction p = po. I emphasise that these functions have 
been arbitrarily chosen for illustrative purposes, though 
with care to make sure that they have the appropriate 
limiting behaviours. In Eq. (|4]) a and /3 are 'stiffness' 
coefficients. They are related to derivatives of the func- 
tions at the fixed point by a = — nA'/A and [3 ^ q' . For 
results presented below I use a — 2 and /3 = 10 which 
allows for stochastic fluctuations of moderate amplitude, 
balanced between the two cell types. I have checked that 
the results are qualitatively insensitive to this choice. 

Fig. [DJb) shows the phase space portrait and a repre- 
sentative stochastic trajectory for this ASPC model with 
a fixed point chosen to lie at po = 0.22 and a target popu- 
lation size TVq = Ano = 200. Again A does not need to be 
explicitly specified though in this case it could be inter- 
preted as representing the area of influence of diffusible 
intercellular signalling factors. In contrast to Fig. [Ha), 
there is an isolated stable fixed point, whose basin of at- 
traction extends to cover the whole plane. The stochastic 



trajectory is strongly localised to the vicinity of the fixed 
point. It is clear that this behaviour should be generic 
to a wide class of ASPC models since the existence of an 
isolated stable fixed point is a structurally stable feature 
of the dynamics. The fixed point in a model of this type 
represents homeostasis in several ways. Firstly, if the 
cell populations are perturbed, they will tend to return 
to the original fixed point. Secondly, fluctuations in the 
cell populations will be limited to the vicinity of the fixed 
point. Thirdly, the model itself can be perturbed, for ex- 
ample by the inclusion of stem cells, without leading to 
a qualitative change in behaviour. 

A key requirement of such models is that they re- 
tain compatibility with the keratinocyte labelling data 
of Clayton et al. 0, Q . An (admittedly mean-field) argu- 
ment that this is true can be made as follows. Imagine 
labelling a small fraction of keratinocytes. If the label is 
passive, proliferation of labelled cells will be determined 
by fixed parameter values, corresponding to the homeo- 
static fixed point. In particular the symmetric division 
pathways of labelled progenitor cells will be automati- 
cally balanced. 

I now turn to the implications for cancer modelling. 
As outlined in the introduction, the ASPC models might 
offer some interesting insights into carcinogenesis, which 
can be considered as unlimited cell proliferation caused 
either by a loss of stability of the homeostatic fixed 
point, or by a large but rare stochastic fiuctuation caus- 
ing the system to exit the homeostatic basin of attrac- 
tion. Whilst the idea that cancer arises from an instabil- 
ity in the underlying cell population dynamics is rather 
old [l^, [l3| , the second possibility is very intriguing and 
has apparently not been hitherto considered. My studies 
of models comprising A, A* and B cells, with a pro- 
cess A ^ A* representing a somatic mutation (c. /. Klein 
et al. Q), shows that the phenomenon of 'homeostatic 
metastability' can easily be observed. However the re- 
sulting three-dimensional phase space portraits are tricky 
to represent and analyse. To illustrate the mechanism of 
homeostatic metastability therefore, I return to the origi- 
nal ASPC model, but introduce a re-entrant bias control 
function g(p). Such a model is a prototypical example of 
a system which is near a saddle-node bifurcation. 

An example of a re-entrant q{p) is given by inserting 
an extra factor (pi — p)/{pi — po) in the argument to the 
tanh function in Eq. (|4|) . This model has a stable node at 
p = Po and a saddle dX p ~ pi. The saddle-node bifurca- 
tion is approached as Ap ~ pi — po vanishes. The phase 
space portrait and a representative stochastic trajectory 
for this ASPC model are shown in Fig.[IJc) for Ap — 0.04 
(other parameters as for the original ASPC model). The 
homeostatic basin of attraction is now confined to the 
lower left region. The saddle lies on the homeostatic 
basin boundary. The simulations show that the system 
inevitably escapes from the homeostatic basin, typically 
in the vicinity of the saddle. After this the cell popula- 




FIG. 1: (color online). Phase space portraits for (a) the SPC 
model, (b) an autoregulating SPC (ASPC) model, and (c) an 
ASPC model exhibiting homeostatic metastability. The axes 
are the progenitor cell fraction, p, and the ratio between the 
current and initial total cell densities, n/no. Thin black lines 
are phase space flows, with the direction indicated by the ar- 
rows. In (a) the thick blue line is the line of fixed points of the 
SPC model. In (b) and (c) the filled blue circles are homeo- 
static stable fixed points (nodes). In (c) the open blue circle 
is an unstable fixed point (a saddle), lying on the homeostatic 
basin boundary shown as a dashed blue line. Jagged red lines 
are representative stochastic trajectories. 



FIG. 2: (color online). Data collapse of homeostatic escape 
rate u as a function of NoAp^. Data were collected for various 
iVo and (colored) for Ap = 0.030(0.005)0.055. 



tions grow without limit. 

To characterise the escape event, I generate a large set 
of trajectories and compute an escape rate u from the 
probability of remaining in the homeostatic basin [l8| . 
Fig. [5] shows that u decreases approximately exponen- 
tially with NqAp^ [l^. The sensitive dependence on the 
target population size No is reminiscent of system size 
effects found for other non-equilibrium phase transitions 
[iol . The sensitivity to the distance Ap from the saddle- 
node bifurcation may reflect a quadratic dependence of 
the height of an 'action' barrier [2lj . 

The concept of homeostatic metastability fits neatly 
into multi-stage models which are widely used to inter- 
pret cancer epidemiology 0, @] ■ In a multi-stage model, 
cell lineages slowly transit through one or more pre- 
cancerous stages before entering a final cancerous stage. 
The idea presented here is that homeostatic escape could 
be the final slow step after one or more somatic muta- 
tions have taken place, bringing, for instance, the system 
close to a saddle- node bifurcation. Note that, extrapo- 
lating from Fig. [21 the homeostatic escape rate can easily 
be comparable to multi-stage transition rates which are 
of the order 10~^-10~^ yr~^ In this context it is 

intriguing to note that the initial step in the develop- 
ment of skin cancer often appears to involve a mutation 
in the ras signalling pathway [23] , which is important for 
controlling cell proliferation. 

The mechanism of epithelial renewal in the lungs has 
much in common with that of the skin, although the 
turnover time may be somewhat longer 2J]. Lung can- 
cer might therefore be expected to have many features 
in common with skin cancer. A long-standing puzzle in 
multi-stage models of lung cancer has been the apparent 
indifference of the rate of the final step to the presence 



of a carcinogen (tobacco smoke) [2^, 2^ 27 1. This epi- 
demiology has been interpreted as indicating that a non- 
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mutagenic mechanism might be at work. One s ugg estion 
is that the final step is epigenetic in nature [26|. An- 
other suggestion is that the final step somehow involves 
signalling . Homeostatic metastability can be seen as 
a specific example of the latter signalling mechanism. As 
Frank challenges though [sj, any such interpretation of 
the epidemiology should be supported by other evidence. 
Experiments that directly test the mechanism of homeo- 
static escape are therefore very desirable! A central fea- 
ture of the present model, around which experiments may 
perhaps be designed, is that the cells themselves do not 
undergo any change if the system escapes from home- 
ostasis; only the micro-environment changes. 

Some general points can be made about directions for 
future research. Firstly, many of the arguments presented 
here are mean-field in nature. This key point was well 
appreciated by Klein et al. Q, who showed that many 
mean-field results still hold in a two-dimensional version 
of the original SPC model. The extension of the present 
ASPC models to fully-fledged two dimensional models is 
clearly a crucial next step. Another direction in which 
progress could be made is to improve the representation 
of the biology, for exarnple moving to multi-scale [1^, 
or agent-based models [29|, which capture the detailed 
biology of the intercellular signals, and also the essential 
stochastic nature of individual cell fates. 

A generic conclusion of the present study is that it 
may not be valid to examine just the deterministic conse- 
quences of somatic mutations, since rare stochastic events 
may occur at comparable rates. This makes the task 
of examining the behaviour of more biologically detailed 
models rather formidable. Brute force methods {i. e. lots 
of very long simulations) have been used in the present 
paper since the underlying stochastic processes are rather 
simple. This may not be possible for more complex mod- 
els. Instead, it may be necessary to bring to bear more 
sophisticated techniques such as transition path sampling 
[30| , or forward-flux sampling 31[ . 

I thank Rosalind Allen, Mike Gates and Martin Evans 
for helpful discussions and encouragement. 
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